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FICTITIOUS  DOMAIN  METHODS 
FOR  VISCOUS  FLOW  SIMULATION 

Roland  Glowinski1,  Anthony  J.  Kearsley2, 
Tsorng-Whay  Pan3,  Jacques  Periaux4 


Abstract 

We  discuss  the  fictitious  domain  solution  of  the  Navier-Stokes  equations 
modelling  unsteady  incompressible  viscous  flow.  The  method  is  based  on 
a  Lagrange  multiplier  treatment  of  the  boundary  conditions  to  be  satis¬ 
fied  and  is  particularly  well  suited  to  the  treatment  of  no-slip  boundary 
conditions.  This  approach  allows  the  use  of  structured  meshes  and  fast 
specialized  solvers  for  problems  on  complicated  geometries.  Another  inter¬ 
esting  feature  of  the  fictitious  domain  approach  is  that  it  allows  the  solution 
of  optimal  shape  problems  without  regriding.  The  resulting  methodology 
is  applied  to  the  solution  of  flow  problems  including  external  viscous  flow 
past  oscillating  rigid  body  and  vortex  dynamics  of  two-dimensional  flow 
modelled  by  the  incompressible  Navier-Stokes  equations  and  then  to  an 
optimal  shape  problem  for  Stokes  and  Navier-Stokes  flows. 


1.  INTRODUCTION 

Fictitious  domain  methods  for  partial  differential  equations  are  showing  interest¬ 
ing  possibilities  for  solving  complicated  problems  motivated  by  applications  from 
science  and  engineering  (see,  for  example,  [1]  and  [2]  for  some  impressive  illustra¬ 
tions  of  the  above  statement).  The  main  reason  for  the  popularity  of  fictitious 
domain  methods  (sometimes  called  domain  embedding  methods ;  cf.  [3])  is  that 
they  allow  the  use  of  fairly  structured  meshes  on  a  simple  shape  auxiliary  domain 
containing  the  actual  one,  therefore  allowing  the  use  of  fast  solvers. 

In  this  article  which  follows  [4-6],  we  consider  the  fictitious  domain  solution 
of  the  Navier-Stokes  equations  modelling  the  unsteady  incompressible  Newtonian 
viscous  fluids  and  apply  the  resulting  methodology  to  the  solution  of  optimal  shape 
problems  for  Stokes  and  Navier-Stokes  flows. 

The  principle  of  fictitious  domain  methods  discussed  here  is  to  solve  the  problem 
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in  a  large  domain  (containing  the  actual  one)  with  a  very  simple  shape,  the  fictitious 
domain,  and  to  impose  the  boundary  condition  by  the  introduction  of  a  Lagrange 
multiplier  on  the  actual  boundary.  Its  advantage  is  to  allow  the  numerical  treatment 
on  a  fixed  structured  mesh,  independent  of  the  actual  boundary  of  bodies  which 
might  be  moving.  Thus  the  time-consuming  construction  of  a  boundary-fitted  mesh 
for  each  different  position  of  moving  bodies  can  be  skipped. 

The  methods  discussed  here  go  far  beyond  the  related  work  in  [7]  where  only  the 
steady  Stokes  problem  was  considered  (in  the  particular  case  where  the  boundary 
of  the  actual  domain  is  compatible  with  the  finite  element  mesh  used  in  the  aux¬ 
iliary  domain).  The  methods  described  in  the  following  sections  do  not  require  a 
strong  coupling  between  the  actual  boundary  discretization  and  the  grid  used  in  the 
auxiliary  domain.  It  also  relies  on  the  splitting  methods  described  in,  e.g.,  [8] — [12]; 
with  these  methods  one  can  decouple  the  numerical  treatments  of  the  incompress¬ 
ibility  and  of  the  advection,  and  take  advantage  of  this  fact  to  use  the  embedding 
approach  in  the  (linear)  incompressibility  step  only,  the  advection  being  treated  in 
the  larger  domain  without  concern  -in  some  sense-  for  the  actual  boundary. 

The  content  of  this  article  is  as  follows:  In  Section  2  we  discuss  fictitious  domain 
methods  for  Dirichlet  problems  and  its  numerical  solution  by  combination  of  do¬ 
main  decomposition  and  fictitious  domain  methods;  then  in  Section  3  we  consider 
fictitious  domain  methods  for  incompressible  Navier-Stokes  equations.  In  Section 
4  the  simulation  of  external  incompressible  viscous  flow  past  an  oscillating  rigid 
body  modelled  by  Navier-Stokes  equations  is  discussed.  In  Section  5  we  apply  the 
methods  of  the  above  sections  to  study  the  vortex  dynamics  of  a  two-dimensional 
viscous  flow  taking  place  in  a  disk.  In  Section  6  we  address  the  fictitious  domain 
solution  of  optimal  shape  problems  for  Stokes  and  Navier-Stokes  flows.  Finally,  in 
Section  7  we  conclude  the  paper  with  some  observations  and  comments  on  future 
work. 

2.  FICTITIOUS  DOMAIN  METHODS  FOR  THE  DIRICHLET 
PROBLEM 

2.1  A  model  problem 


Fig.  1. 


We  consider  the  following  Dirichlet  problem: 


au  —  uAu  —  f  in 

(1) 

u  =  g0  on  7, 

(2) 

«  =  Sion  r, 

(3) 

where  Q  is  a  “box”  shaped  domain  in  R d(d>  1),  u>  is  a  bounded  domain  in  Rd  ( d  >  1) 
such  that  to  is  contained  in  fl  (e.g.,  see  Fig.  1),  T  (resp.,  7)  is  the  boundary  dfi 
(resp.,  du>),  we  suppose  that  7  is  smooth;  a  >  0  and  v  >  0;  finally,  /,  g0 ,  and  gx 
are  given  functions  defined  over  7  and  T,  respectively.  If  /,  g0,  and  gx  are 

smooth  enough,  problem  (1)  (3)  has  a  unique  solution. 

A  fictitious  domain  method  was  already  proposed  for  problem  (l)-(3)  in  [4],  For 
simplicity,  we  shall  assume  that  f  E  L2(fl\u).  We  imbed  fl\u>  in  ft  and  define 

Vgi  =  {u|v  €  Hx( 0),u  =  gx  on  T}.  (4) 


We  then  have  equivalence  between  (l)  (3)  and  the  following  saddle-point  problem: 
Find  {u,  A}  €  Vgi  x  L2( 7)  such  that 


/  (auv  +  vVu  •  Vu)  dx  =  /  fv  dx  +  /  Xv  dy,  \/v  E  Hq  (fl), 

v  n  i/fi  J'y 

/  -  9o )  d'y  =  0 ,  V/16  L?{ 7), 


(5) 

(6) 


where  /  is  a  L2(Sl)-extension  of  /  and  satisfies  =  /.  It  can  be  shown  that 

u|^-  is  the  solution  of  problem  (l)-(3)  and  «  =  j  on  7.  In  (5)-(6),  more  precisely 
we  have  that  the  multiplier  A  is  equal  to  i/[dtt/dn]|7  (i.e.,  u  times  the  jump  of  the 
normal  derivative  of  u  at  7). 

Problem  (5)-(6)  can  be  solved  by  the  general  conjugate  gradient  methodology 
proposed  in  [4].  For  the  cases  where  u  C  R2  with  a  smooth  boundary  7,  we  have 
obtained  by  Fourier  Analysis  a  quasi-optimal  preconditioned  conjugate  gradient 
algorithm  for  the  solution  of  (5)-(6)  (see  [4]). 

2.2  Domain  decomposition/fictitious  domain  approach 

Motivated  by  computation  with  nonmatching  grids  on  MIMD  machines,  we  con¬ 
sider  the  numerical  solution  of  Dirichlet  problems  by  a  combination  of  domain 
decomposition  and  fictitious  domain  methods.  We  take  advantage  of  the  fact  that 
the  Steklov-Poincare  operators  associated  with  the  subdomain  interfaces  and  with 
the  fictitious  domain  treatment  of  internal  boundaries  have  very  similar  properties. 
We  use  these  properties  to  derive  fast  solution  methods  of  the  conjugate  gradient 
type  with  good  parallelization  properties  which  simultaneously  force  the  matching 
at  subdomain  interfaces  and  the  actual  boundary  conditions. 


Fig.  2. 


We  consider  the  case  where  u>  is  a  bounded  domain,  and  a  three  subdomain 
decomposition  like  the  one  in  Fig.  2  where  ft  =  fix  U  f22  U  ^3;  we  denote  by 
7i2  (resp.,  723)  the  interface  between  Oi  and  fi2  (resp.,  f22  and  f23),  and  by  7 
the  boundary  of  w;  and  let  T;  =  T  D  dfli  for  i  =  1,2,  3.  To  obtain  the  equivalent 
fictitious  formulation  of  (l)-(3),  we  imbed  fl\d)  in  fl  and  define  the  following  spaces 

Vgi  =  iv\v  e  FTx( fli ),v  =  gx  on  Tj,  Vq  =  {t>|u  G  H1{ fli),v  =  0  on  T;}, 
for  i  =  1,  2,  3  and  let 

v  =  Vg\  x  V2  x  V* ,  V„  =  V0'  x  V02  x  Vi,  A  =  L2(7)  x  L2(7l2)  x  £2(723). 
Problem  (l)-(3)  is  equivalent  to  the  following  system 


Find  (ux,U2,  u3)  G  V ,  (A,  Ai2,  A23)  G  A  such  that 


/  {auiVi  +  vVui  ■  Vvi)  dx  =  /  fvidx+  /  A^t^  —  17)  ds  (7) 

j_l  J£ii  d 

+  /  A23(v3  -  u2)  ds  +  /  Av2d7,  V(t»i,t;2,v3)  G  V0, 

7 -/no  J  'Y 


/  -  <7o)  dy  =  0,  VgG  L2(l), 

h 

(8) 

/  Mi2(«2  -  «i)d7  =  0,  V/*i2  G  L2(7i2), 

(9) 

'712 

/  M23(«3  -  M2)  «*7  =  0,  V/u23  G  L2( 723), 

J-m 

(10) 

where  /  is  a  L2(fl)  extension  of  /.  We  have  equivalence  in  the  sense  that  if  relations 
(?)  (10)  hold  then  Ui  —  for  i  —  1,2,3,  where  u  is  solution  of  (5)-(6),  and 
conversely. 

Due  to  the  combination  of  the  two  methods,  there  are  two  Lagrange  multipliers 
in  (7)— (10) .  The  function  A  which  is  a  Lagrange  multiplier  associated  with  the 
boundary  condition  u  =  g0  on  7  is  essentially  the  jump  of  vidujdn )  at  7,  and  the 
function  A12  (resp.,  723)  which  can  be  viewed  as  a  Lagrange  multiplier  associated 
with  the  interface  boundary  condition  Ui  =  u2  (resp.,  u2  =  u3)  on  712  (resp.,  723)  is 
nothing  but  the  function  u(du/dm)\112  =  -u(du/dn2) |7la  (resp.,  v(du/dn2)\y23  = 
-p(du/dn3)\l23)  where  nl  is  the  outward  normal  unit  vector  of  O,  at  interfaces  712 
and  723. 

There  are  three  different  conjugate  gradient  algorithms  which  have  different  par¬ 
allelization  properties  and  can  be  parallelized  on  MIMD  machines  [13].  Here  we 
apply  the  one  called  the  one  shot  method  to  solve  the  saddle— point  system  (7)— (10). 
In  [14],  the  numerical  experiments  of  the  one  shot  method  to  two-dimensional  and 
three-dimensional  test  problems  have  been  done  on  a  KSR1  of  MIMD  type. 

The  one  shot  method  is  the  following: 


A0  =  (A0,  A?2,  A23)  €  A  given ; 


solve 

Find  (u°,  «2,  W3)  €  V,  such  that 

^2  f  (au°i  Vi  +  uVui  •  Vvi)  dx  =  l  fvidx+  f  \°12(v2  -  ds 
i  =  l  “  £li  j  — 1  J n, 

+  /  Aj 3(v3-v2)ds  +  /  X°v2dj,  \/(vi,v2,v3)  e  Vo, 

d  723  J  7 

sat  g°  =  (g°,  g®2,  g®3)  —  ((u°  -  #o)|7,  (u2  —  «?)|7l2,  (“3  -  ^2) I-X23 )? 

and  set  w°  =  {w°,  w°12,  w°3)  =  (g°,  g°2,  g°3). 

For  n>  0,  knowing  Xn,  gn,  w",  compute  An+1,  gn+1,  wn+1  as  follows: 
solve 

Find  (ui,u2,u3)  G  Vo,  such  that 

Y.  /  (aufvi  +  uVu^  ■  Vvi)dx  =  j  w^2(v2  -  Vi)  ds 

i= 1  d^i  J  712 

+  /  w23(v3  -  v2)  ds  +  /  wnv2  d'y ,  V(ui,  v2,  v3)  €  Vo, 

d 723  J  7 

and  set  g"  =  (gn,g^2,g^3)  =  («£ I7,  K  -  «i*)l7.2>  («a  - 
We  then  compute 

Jj  \dn\2  d'y  +  fji2  \g?2\2  ds  +  fj23  \glf3\2  ds 


Pn  — 


J7  gV  dj  +  fyiy  g«2w?2  ds  +  /  gf3wf3  ds  ’ 


(11) 


(12) 


(13) 


(14) 


and  set 


(15) 

(16) 
(17) 


An+1  =  An  _ 

u?+1  =  u7  -  Pnrf},  for  i  =  1,  2, 3, 

gn+i  =  gn  _  pngn_ 


If  ||gn+1||/||g°||  <  e,  take  A  =  An+1,  iti  =  «i+1,  u2  =  «£+\  anc^  M3  =  «3+1/  */ 

not  compute 

T„  =  l|g"+1||2/l!gn||2,  (18) 

and  set  wn+1  =  gn+1  +  7nwn. 

Do  n  =  n  +  1  and  go  to  (13). 

In  algorithm  (11)-(18),  ||g||2  =  /7|p|2^7  +  /7l2  |0i2|2ds  +  /723  |.92s|2  for  any 
g  =  (9,912,923)  £  A.  We  have  taken  e  =  10-8  in  the  stopping  test.  When  u  C  R2 
with  smooth  boundary,  a  preconditioned  one  shot  method  has  been  discussed  in 

[13]- 

2.3  Numerical  experiments 

We  consider  the  following  test  problem: 

—  Ait  =  0  in  (19) 

it  =  0  on  7,  (20) 

it  =  — tccos#  +  ysin#  on  T,  (21) 

where  fl  =  (—0.5,  0.5)  x  (—0.5, 0.5),  to  is  a  NACA0012  airfoil  centered  at  (0,0)  with 
zero  degree  of  angle  of  attack  whose  chord  length  is  0.35,  and  9  is  an  incident  angle. 
In  numerical  experiments,  U  U  ^3  where  fix  =  (—0.5, 0.5)  x  (0.125,  0.5), 

tt2  =  (-0.5,  0.5)  x  (-0.125,  0.125),  =  (-0.5,  0.5)  x  (-0.5,  -0.125). 

The  finite  dimensional  spaces  V*ih  and  V^h  which  approximate  Vxgx  and  Vq  ,  re¬ 
spectively,  for  i  =  1,2,3  are  defined  as  follows: 

Vgih  =  {vh\vh  €  C°m,vh  -  gih  on  Tz,vh\T  e  Pi,  VT  €  Iff), 
voh  =  ivh\vh  e  C°(S),i),vh  =  0  on  Vi,vh\T  €  Pi,  VT  e  Tff}, 

where  gm  is  an  approximation  of  gi  given  by  (21),  is  a  triangulations  (see  Fig. 
3)  of  Qi  for  i  =  1,  2,3  and  Pi  is  the  space  of  the  polynomials  in  xi,x2  of  degree 
<  1.  For  the  above  finite  dimensional  spaces,  we  can  use  different  mesh  size  in  each 
subdomain  (see  Fig.  3). 


The  spaces  L2( 7),  L2( 712)  and  L2( 723)  are  approximated  by  the  finite  dimen¬ 
sional  space  Aft,  Afti2,  and  Aft23,  respectively,  as  follows: 

A ft  =  {Mft|/ift  €  L°°( 7),  ph  is  constant  on  the  segment  joining 
2  consecutive  mesh  points  on  7}, 

Afti2  =  {ph\hh  €  L°°( 7i2),  ph  is  constant  on  the  segment  joining 
2  consecutive  mesh  points  on  7i2}, 

Aft23  =  {ph.\ph  €  T°°( 723),  /^ft  is  constant  on  the  segment  joining 
2  consecutive  mesh  points  on  723}. 

A  particular  choice  of  mesh  points  on  7  is  visualized  in  Fig.  3.  The  mesh  points 
on  7i2  and  723  are  the  midpoints  of  the  edges  located  on  the  interfaces  of  the  finest 
triaugulation  [15]  (e.g.,  see  Fig.  4). 


Fig.  3:  Example  of  meshes  on  fix,  ft2,  and  fl3  with  mesh  sizes  h  =  1/8, 
h  =  1/32,  and  h  =  1/8  respectively.  Mesh  points  on  7  are  marked  by  “*”. 


Fig.  4:  Example  of  mesh  on  interface  between  subdomain  with  nonmatch¬ 
ing  grids.  Mesh  points  on  interface  are  marked  by  “1” 


The  solution  of  (19)-(21)  is  the  stream  function  of  the  potential  flow  of  incident 
angle  6.  In  the  numerical  experiments,  we  have  taken  9  =  45°  and  the  mesh  sizes 
are  hi  =  1/128,  h2  =  1/512,  and  h%  —  1/128  for  the  subdomains  S7i,  Cl2,  and  O3 
respectively.  We  applied  the  one  shot  method  with  nonmatching  grids  to  (19)  (21) 
and  the  corresponding  numerical  results  are  shown  in  Fig.  5.  Those  obtained  by  the 
fictitious  domain  method  with  uniform  grid  of  mesh  size  h  =  1/512  are  also  shown 
in  Fig.  5.  With  the  one  shot  method  and  nonmatching  grids,  we  can  successfully 
use  fine  grid  in  the  region  close  to  the  body  and  coarse  grid  in  the  region  far  from 
the  body.  We  observe  the  excellent  agreement  between  these  results. 


Fig.  5:  The  left  figure  shows  the  streamlines  obtained  by  the  one  shot 
method  with  nonmatching  grids  and  h1  =  1/128,  h2  =  1/512,  and  h3  = 
1/128  for  the  subdomains  Hi,  122 ,  and  O3  shown  in  Fig.  3.  The  right 
figure  shows  the  streamlines  obtained  by  the  fictitious  domain  method  with 
uniform  grid  of  mesh  size  h  =  1/512. 

3.  FICTITIOUS  DOMAIN  METHODS  FOR  THE  NAVIER-STOKES 
EQUATIONS 

3.1  The  Navier-Stokes  Equations 


r„ 


Fig.  6. 


In  [5,  6],  we  have  considered  external  incompressible  viscous  flow  modeled  by 
the  Navier— Stokes  equations  with  either  Dirichlet  boundary  conditions  or  Neumann 
boundary  conditions  downstream.  Using  the  notation  of  Fig.  6,  we  consider  the 
following  problem 


—  -i/Au  +  (u'  V)u  +  Vp  =  f  intl\u(t),  (22) 

V  •  u  =  0  in  12  \  u>(t),  (23) 

u(x,0)  =  u0(x),  x  €  (with  V  •  u0  =  0),  (24) 

dti 

u  =  g0(t)  onF0,  u—-np  =  g1(t)onTi,  (25) 

u  =  g2(t)  ony(t).  (26) 


In  (22)-(26),  tt  is  bounded  domains  in  I d(d  >  2),  co(t)  could  be  a  moving  rigid 
body  in  Rd(d  >  2)  (see  Fig.  6),  T  (resp.,  y(t))  is  the  boundary  of  Q  (resp.,  ui(t)) 
with  r  =  r0  UTi,  r0nri  =  0,  and  Jr  dT  >  0,  n  is  the  outer  normal  unit  vector 
at  Ti,  u  =  {uiYjZf  is  the  flow  velocity,  p  is  the  pressure,  f  is  a  density  of  external 


dwi 


forces,  v{>  0)  is  a  viscosity  parameter,  and  (v  •  V)w  =  Vj - -} \~f. 

^  dxj 

To  obtain  the  equivalent  fictitious  domain  formulation  for  the  Navier-Stokes 
equations,  we  embed  fi  \  u>(t)  in  and  define 


Vgo  =  (vlv  e  (H1(n))d,  v  =  g0  on  r0},  v0  =  {v|v  e  (h1^))4*,  v  =  oon  r0}, 

and  A (t)  =  (L2(7(t)))d.  We  observe  that  if  U0  is  an  extension  of  u0  with  V-Uo  =  0 
in  Q,  and  if  f  is  an  extension  of  f,  we  have  equivalence  between  (22)-(26)  and  the 
following  problem: 


For  t  >  0,  find  U(i)  e  Vg0,  P(t)  e  L2(J7),  A (t)  6  A (t)  such  that 

f  d  U  ,  r 

/  — —  •  v  dx  +  v  /  VU  •  V v  dx 

Ju  Ju 

+  [  (U-V)U-vdx—  f  PV  •  v  dx 

Jn  Jn 

=  /  f  •  v  dx  +  /  g!  •  v  dT  +  /  A  •  V  dj,  Vv  e  V0,  a.e.  t  >  0, 

J  a  dri 


(27) 


V  •  U(f)  =  0  in  ft, 

(28) 

U(x,  0)  =  U0(x),  x  €  ft,  (with  V  ■  U0  =  0), 

(29) 

U(t)  =  g 2(f)  on  7(f), 

(30) 

in  the  sense  that  =  u,  P\n\ZJ(t)  =  P'  Above,  we  have  used  the  notation 

4>{t)  for  the  function  x  — >  d>(x,  t). 

Concerning  the  multiplier  A,  its  interpretation  is  simple  since  it  is  equal  to  the 
jump  of  v(d\J /dn)  -  nP  at  7 (t).  A  closely  related  approach  (limited  to  the  steady 


Stokes  problem)  is  discussed  in  [7].  We  observe  that  the  effect  of  the  actual  geometry 
is  concentrated  on  J A  •  v  dj  in  the  right-hand-side  of  (27),  and  on  (30). 

3.2  Time  discretization  by  operator  splitting 


To  solve  (27)-(30),  we  shall  consider  a  time  discretization  by  an  operator  splitting 
method ,  like  the  ones  discussed  in,  e.g.,  [8]— [12] .  With  these  methods  we  are  able  to 
decouple  the  nonlinearity  and  the  incompressibility  in  the  Navier-Stokes/fictitious 
domain  problems  (27)-(30).  In  the  following,  we  apply  the  9-scheme  (cf.  [12]) 
to  (27)-(30)  with  time  step  At  >  0.  Let  As  =  (L2(7(sAf)))d  and  denote  by  <j)s 
either  an  approximation  of  or  the  function  cf>(sAt).  We  obtain  the  following  time 
discretization  scheme: 

U°  =  Uo  is  given ;  (31) 

for  n  >  0,  knowing  Un,  find.  Un+(9  G  Vg«+# ,  Pn+9  g  L2(Q),  An+e  G  An+6>  such  that 


t 

■If, 


Un+e  -  Ur 


n  9  At 

=  f  f  n+*-vdx-  f  (Un  •  V)Un  •  v  dx  +  [  A  n+e-vd7 

Jo,  J  n 

-pu  I  VUn  •  Vv  dx  +  f  (a^+6  +  /3g")  •  v  dT,  Vv  G  V0, 
■Isa  J  rt 


•  v  dx  T  olv  j  VUn+e-Vvdx-  f  Pn+eVvdx 

J  £1  i/ 


(32) 


V  •  Un+<?  =  0  in  O,  (33) 

Un+*  =  g£+e  on  7n+e;  (34) 

next,  find  Un+1-e  G  V_n+i-«  such  that 

So 


[  L 


pri+l-e  _  jjn+0 


,  .  •  vdx  +  /?i/  /  VUn+1"e  •  Vvdx 

(1  -  29)  At  H  Ja 

+  f  (uri+1“e  •  V)Un+1-*  •  V  dx 
J  n 


=  [  f n+1~9-vdx+  f  Xn+e  •  v d7  +  f  Prt+*V-vdx 
Jsa  J^n+e  Jn 

—olv  [  VUn+e  •  Vv  dx  +  f  (a^+e  +  0f£+1-e)  •  v  dT,  Vv  G  V0; 
^  Jo  JTi 

finally,  find  Un+1  G  V  ™+i ,  Pn+1  G  L2(tt),  A"+1  G  An+1  such  that 

So 


(35) 


f  / 

J  o 


■jjn+l  _  pjn+1- 


9  At 


j 

•In 


•vdx  +  cw/  /  VUn+1  •  Vv  dx  —  /  Pn+1V-vdx 


/■ 


[  f n+1-vdx-  [  (U"+1-®  ■  V)Un+1"e  ■  vdx+  [  An+1  •  v  d7  (36) 
dn  dn  i7»+i 


f  VUn 
d  o 


+1“e  •  Vvdx  + 


/  («grx 

dr! 


+  /3gi+1_0)-vdr,  VvG  Vo, 


v 


V  •  Un+1  =  0  in  ft,  (37) 

Un+1  =  S2+1  on  7n+1,  (38) 

where  a  +  (3  =  1,0  <  a,  )3  <  1  and  0  <  6  <  1/2.  With  the  choice  of  0  =  1  -  l/y/2, 
a  —  2  —  \[2  and  f3  =  y/2  —  1,  the  time  discretization  seems  to  be  unconditionally 
stable  (see  [12]). 

In  Section  3.3  the  conjugate  gradient  solution  of  the  Stokes/ fictitious  domain 
problems  (32)-(34)  and  (36)-(38)  will  be  discussed.  Concerning  problem  (35)  it  is 
worth  noticing  that  we  have  been  taking  advantage  of  the  time  discretization  by 
operator  splitting  to  treat  the  advection  in  the  larger  domain  ft  without  concern 
-in  some  sense-  for  the  constraint  u  =  g  at  7.  Problem  (35)  can  be  solved  by  least- 
squares  methods  [12]  and  is  also  well  suited  to  solution  methods  based  on  higher 
order  upwinding  on  regular  meshes,  or  on  the  backward  method  of  characteristics 
(see,  e.g.,  [16]). 

3.3  Iterative  solution  of  the  Stokes/fictitious  domain  problem 

Problems  (32)-(34)  and  (36)-(38)  are  particular  cases  of  the  following  Stokes/fic¬ 
titious  domain  problem: 


(  Find  U  G  Vg0,  P  G  L2(ft),  A  G  A  such  that 


/  U  •  v  dx  +  v  / 

«/  J  £7 


dyi  +  v  /  VU-Vvdx 


=  /  f  •  v  dx  +  /  A-vd7  + 
^  in  J-y 


/  PV  •  v  dx 
in 

[  gi  •  vdr,  Vv  6  v0, 

J  Fi 


V  •  U  =  0  in  ft, 
U  =  g 2  on  7, 


(39) 

(40) 

(41) 


where,  in  (39),  a(>  0)  is  the  reciprocal  of  a  partial  time  step.  In  (39)-(41),  P 
(resp.,  A)  appears  to  be  a  Lagrange  multiplier  associated  with  (40)  (resp.,  (41)). 

We  can  solve  the  above  saddle  point  system  (39)-(4l)  by  a  conjugate  gradient 
algorithm  (called  the  one  shot  method)  driven  by  the  pressure  P  and  the  multiplier 
A,  simultaneously.  Let  us  consider  a  bilinear  form  6(-,  •),  symmetric  and  strongly 
elliptic  over  A  (we  may  choose  b( A,p)  =  /  A  ■  fi  d-y,  VA,  //  e  A).  The  following 
algorithm  is  the  one  shot  method  driven  by  the  multipliers  P  and  A: 

{P°,  A0}  e  L2(ft)  x  A  given ;  (42) 

solve  the  following  Dirichlet  problem 

{  Find  U°  G  Vg0 ,  such  that 

a  [  U°  •  v  dx  +  u  [  VU°  •  Vv  dx  =  f  f  •  v  dx+ 

Jn  in  ./  n 

/  A0  •  v  d7  +  f  P°  V  •  v  dx  -f  f  gi  •  v  dT,  Vv  G  V0, 
i  7  in  Jty 


< 


(43) 


set  r°  =  V  •  U°,  =  (U°  —  g2)|7  and  define  g°  =  {5°,  g^}  as  follows 


0  iO  ,  0 

9i  =  a(P  +  vr  1, 


with  (f>°  the  solution  of 


—A (f>°  =  in  fi, 
d(f>° 


dn 


0  on  To;  4>°  =  0  on  IA; 


Find  g2  €  A  such  that 
&(g2»  /*)  =  /  r2  •  M d7,  V/t  €  A. 

1/  y 

We  take  w°  =  {«,?,  w°}  =  g°}. 

Then  for  n  >  0,  assuming  that  { Pn ,  An},  Un,  r",  ,  wn,  gn  are  known, 

{Pn+1,  An+1},  Un+\  r”+1,  r£+l,  wn+1,  gn+1  as  follows: 


solve 


(  Find  Uri  €  Vq,  such  that 


J 

Jn 


a  I  Un-vdx  +  v  /  VUn-Vv<fx 

Jn 

I  w2  ■  vd'y  +  (  w™  V  •  vdx,  Vv  6  Vo; 

J  7  J  n 

set  f  %  =  V  •  tJn,  fj  =  Un|7  and  define  g"  =  {5”,  g£  }  as  follows: 

g^ar  +  vr?, 

with  f>n  the  solution  of 

—A 4>n  —  frf  in  fl, 

dfi 


dn 


=  0  on  r0 ;  (j)n  =  0  on  T 1] 


Find  g2  G  A  swc/i  t/iai 


^(g2>M)  =  [ 

J  7 


r2  •  V/z  e  A. 


We  then  compute 

ft.  =  (Jn  Ac  +  /7  r?  •  g?  rfy)/(Jn  *=?«>?  d*  +  /,  *2’  w5  <*7), 

and  set 

Xjn+i  =Un_ 

P"+l  =  Pn  —  PnW™,  An+1  =  An  -  pnw%, 


-«+l  —  ~n 


rl  -  Pn»l 


y.n+1 


r 2  -  Pnr'2‘: 


(44) 

(45) 

(46) 

compute 

(47) 

(48) 

(49) 

(50) 

(51) 

(52) 

(53) 

(54) 


</l’  +  1  =  J?  -  P„9?,  g';+1  =  gj  -  P,.g!  (55) 

If  (fa  '■i+ISi+I  dx  +  X, rr 1  ■  g?+1  <*T)/(/n  r?s?  dx  +  fy  r°  •  g°  d7)  <  e,  take  P  = 
Pn+1 ,  U  =  Un+1,  and  A  =  An+1.  //  not,  compute 

■Jn  =  (fa  r"+1<+1  dx  +  fy  p;+1  •  g2n+ 1  dy)/(fa  dx  +  fy  rj  •  gj  i7),  (56) 

and  set  wn+1  =  gn,+1  +  7nwn. 

Do  n  =  n  +  1  and  go  back  to  (47). 


4.  FICTITIOUS  DOMAIN  METHOD  FOR  FLOW  PAST  A  MOVING 
RIGID  BODY 


Fig.  7. 

In  this  section  we  consider  an  incompressible  viscous  flow  modelled  by  the  Navier- 
Stokes  equations  past  a  moving  rigid  body,  a  NACA0012,  by  a  fictitious  domain 
method.  In  the  test  problem,  let  S7  =  (—0.625,0.625)  x  (—0.5,  0.5)  (see  Fig.  7) 
and  w  be  a  NACA0012  airfoil  centered  at  (0,  0)  (half-chord  point)  which  is  forced 
to  pitch  about  its  center  with  the  prescribed  angular  displacement  defined  by 

0(t)  —  90  -  6 a  cos  71 1,  (57) 

where  0o  =  6a  =  10°  and  the  chord  length  is  0.125.  Thus  the  angle  of  attack  is 
between  0°  to  20°.  The  boundary  conditions  are  defined  as  follows: 

11=^^  on  T0,  v -  np  =  O  on  T i .  (58) 

When  constructing  finite  dimensional  subspaces,  we  can  skip  the  time-consuming 
construction  of  a  boundary-fitted  mesh  for  each  different  position  of  the  moving 
body  and  use  a  fixed  uniform  structured  mesh  %  which  is  a  triangulation  of  Q 
(see,  e.g.,  Fig.  8).  For  finite  dimensional  subspaces  approximating  Vg0  (here  let 
go  =  (5o>tfo)  =  (1,0))  and  V0,  we  choose 

vgo h  =  {v^|va  e  V,1  x  vfc2},  VQh  =  {vfclvfc  G  Hlh  x  HqH}, 


(59) 


where 


vh  =  {<f>h\<l>h  e  c°(n),0h|r  e  px,  vt  e  Th,<i>h  =  gl0h  on  ro},for  t  =  1,2,  (60) 

Hq h  -  {4>h\(f>h  e  C°(n)>h|r  €  Pi,  VT  €  Th,4>h  —  0  on  r0}.  (61) 

In  (60)  and  (61)  Pi  is  the  space  of  the  polynomials  in  xi,  x2  of  degree  <  1,  and  <7^ 
(resp.,  g2)h)  is  an  approximation  of  g$  (resp.,  g%).  A  traditional  way  of  approximating 
the  pressure  is  to  take  it  in  the  space 

H\h  =  {4>h\4>h  e  c°(fi),  0fc|T  e  Pi,  vt  e  7k}, 

where  7k  is  a  triangulation  twice  coarser  than  7k  Concerning  the  space  Ash  ap¬ 
proximating  As,  we  define  it  by 

A£.  =  {l^hll^h  €  (L°° (■ys))2 ,  Hh  is  constant  on  the  segment 
joining  2  consecutivemesh  points  on  7s}. 

A  particular  choice  for  mesh  points  on  7s  for  different  angles  of  attack  is  visualized 
on  Fig.  8. 


Fig.  8:  Part  of  the  triangulation  of  Q  with  h  =  1/128  and  mesh  points 
marked  by  on  NACA0012  with  10°  (left)  and  20°  (right)  degrees  of 
angle  of  attack. 


Fig.  9:  Yorticity  distribution  and  stream  lines  of  the  initial  velocity  field. 


Fig.  10:  Vorticity  distribution  and  stream  lines  at  time  t  =5,  5.5,  6,  6.5,  7 
(from  top  to  bottom,  angles  of  attack  are  0°,  10°,  20°,  10°,  and  0°)  during 
one  period  of  the  airfoil  motion  at  Reynolds  number=100. 

Using  the  0-scheme,  we  solve  at  each  time  step  two  Stokes/fictitious  domain 
subproblems  by  the  one  shot  method  discussed  in  Section  3  and  one  advection 


diffusion  subproblem  by  a  least-squares/conjugate  gradient  algorithm  [12].  In  the 
numerical  simulation,  the  Reynolds  number  is  100  (taking  the  chord  of  airfoil  as 
characteristic  length),  the  mesh  size  for  the  velocity  field  is  hv  =  1/128,  the  mesh 
size  for  the  pressure  is  hp  —  1/64,  the  time  step  is  At  =  0.005,  and  e  for  the  stopping 
criterion  in  (42)-(56)  is  10-14.  Fig.  9  shows  the  vorticity  distribution  and  stream 
lines  of  the  initial  condition  for  the  simulation.  It  is  a  steady-state  flow  obtained  at 
zero  degree  of  angle  of  attack.  Fig.  10  shows  a  sequence  of  frames  for  the  vorticity 
distribution  and  stream  lines  during  one  period  of  the  airfoil  motion.  From  top  to 
bottom,  angles  of  attack  are  0°,  10°,  20°,  10°,  and  0°. 

5.  AN  APPLICATION  TO  THE  VORTEX  DYNAMICS  OF  A 
TWO-DIMENSIONAL  FLOW 

Motivated  by  [17,  18],  we  would  like  to  study  the  long-time  evolution  of  the 
vortex  dynamics  of  a  two-dimensional  incompressible  viscous  flow  modelled  by 
the  Navier-Stokes  equations.  In  [19],  selective  decay  and  coherent  vortices  in 
two-dimensional  incompressible  turbulence  were  studied  in  a  square  with  periodic 
boundary  conditions.  In  this  section  we  consider  a  two-dimensional  incompressible 
viscous  flow  in  a  disk. 

In  the  test  problem,  let  Q  =  (—1.5, 1.5)  x  (—1.5, 1.5)  and  w  be  a  disk  of  radius 
1  centered  at  (0,0)  (see  Fig.  11).  Here  we  imbed  w  in  fl.  This  approach  is  slightly 
different  from  the  one  discussed  in  Sections  3  and  4;  but  we  can  use  the  same 
fictitious  domain  techniques.  The  Dirichlet  boundary  conditions  on  7  and  I'  are 
defined  as  follows: 

u  =  0on  r  U  7.  (62) 


Fig.  11. 

As  finite  dimensional  subspaces  of  Vg0  (here  g0  =  (<7o,<7o)  =  (0,0))  and  V0,  we 
choose 

vg0t  =  Wlv*  6  Vi  x  Vi),  V„„  =  {vk|vh  e  X 


where 


(63) 


Vl  =  {<!>h\ct>heCQ{n)^h\TePu  VT  e  T^,  4  =  gl0h  on  r},for  i  =  l,2,  (64) 

Hlh  =  €  C°(n)>fc|r  €  Pi,  VT  €  ThAh  =  0  on  T}.  (65) 

In  (64)  and  (65),  Th  is  a  triangulation  of  fi  (see,  e.g.,  Fig.  12),  Pi  is  the  space  of 
the  polynomials  in  Xi,x2  of  degree  <  1  and  g\h  (resp.,  g1h)  is  approximation  of 
gl  (resp.,  Qq).  A  traditional  way  of  approximating  the  pressure  is  to  take  it  in  the 
space 

H\h  =  {Mh  €  C°(Cl)Ah It  e  Pu  VT  e  7k }, 

where  7k  is  a  triangulation  twice  coarser  than  Th- 

Concerning  the  space  Ah  approximating  A,  we  define  it  by 

A h  —  {l^hl^h,  C  (L°°(7))2,  h h  is  constant  on  the  segment 
joining  2  consecutivemesh  points  on  7}. 

A  particular  choice  for  the  mesh  points  of  7  is  visualized  on  Fig.  12. 

In  Fig.  13  we  show  the  vorticity  and  stream  lines  of  the  initial  velocity  field 
in  the  test  problem  (also  obtained  by  the  algorithm  used  for  solving  the  test 
problem).  In  computing  the  initial  condition  shown  in  Fig.  13,  we  have  chosen 
Uq  =  (sin  7ra:  cos  iry ,  —  cos  nx  sin  Try)  as  initial  velocity,  the  boundary  conditions  be¬ 
ing 

ul  =  e-ctUo  on  T  U  7,  (66) 

where  c  =  5000  in  (66).  The  mesh  size  for  the  velocity  field  is  hv  =  1/120,  the  mesh 
size  for  the  pressure  is  hp  =  1/60,  and  the  time  step  is  At  =  0.001.  The  solution 
obtained  after  1000  time  steps  was  retained  the  initial  condition  for  the  final  test 
problem. 


Fig.  12:  Example  of  the  triangulation  of  with  h  —  1/8  and  mesh  points 
marked  by  on  the  boundary  of  disk. 


In  the  numerical  simulation,  the  Reynolds  number  is  10011  (taking  the  diameter 
of  disk  as  characteristic  length  and  the  maximum  speed  of  the  initial  velocity  field 
as  characteristic  speed),  the  mesh  size  for  the  velocity  field  is  hv  =  1/120,  the  mesh 
size  for  the  pressure  is  hp  =  1/60,  and  the  time  step  is  At  =  0.01.  In  Fig.  14  a 
sequence  of  frames  for  the  vorticity  contours  and  stream  lines  at  t  =  30,  170,  and 
557.5  are  shown.  From  t  =  0  to  about  t  =  120,  vortices  are  stretched,  collide  and 
merge.  By  about  t  =  120,  all  major  vortex-concentrations  have  merged  into  two 
strong  counterrotating  vortices;  there  are  also  two  minor  counterrotating  vortices. 
Then  after  a  long  period  of  time  ( t  =  557.5),  there  are  two  counterrotating  vortices 
surviving.  Let  the  flow  kinetic  energy  be  E  =  |  /^(u i)2  +  (u2)2  dx  and  the  enstrophy 
of  flow  in  the  disk  W  =  f  {du^/dxx  —  du±/dx2\2  dx.  The  evolutions  of  E  and 
W  are  shown  in  Fig.  15.  At  t  =  557.5,  E  has  decreased  by  a  factor  of  6.16  x  10”6 
from  its  initial  value  which  was  0.34731  and  W  has  decreased  to  1.04  x  10” 6  from 
its  initial  value  which  was  49.588. 


Fig.  15:  Time  evolution  of  energy  E  (thin  line)  and  enstropy  W  (thick 
line). 

6.  Shape  Optimization  of  the  Boundary 

Optimal  shape  design  problems  that  are  governed  by  partial  differential  equations 
have  become  a  fruitful  source  of  challenging  and  interesting  free  boundary  problems. 
The  study  of  particular  classes  of  these  problems  has  yielded  solutions  to  fascinating 
applications  (see  for  example  [20-22]).  Recently  researchers  have  tried  to  employ 
the  fictitious  domain  approach  to  problems  of  this  type  arising  in  various  fields  of 
science  and  engineering  (see  for  example  [23-26]).  In  this  section  we  discuss  such 
a  strategy.  A  fictitious  domain  approach  is  employed  to  solve  shape  optimization 
problems  for  both  Stokes  and  Navier-Stokes  equations  in  two  space  dimensions. 

The  shape  of  a  symmetric  airfoil  is  given  by  its  boundary  7.  In  our  test  problems 
we  fix  the  chord  length  of  the  airfoil  hence  keeping  the  location  of  both  the  head 
and  the  tail  of  the  airfoil  fixed.  We  refer  to  [24]  for  more  extensive  numerical  results 
of  our  algorithm  on  test  problems  where  the  location  of  the  head  and  the  tail  of  the 
airfoil  are  design  or  control  variables,  and  also  non-symmetric  shapes  were  used  to 
generate  target  pressures.  With  head  and  tail  fixed,  we  parameterize  the  shape  of 


the  airfoil  by  5  parameters,  say  Kt  where  i  =  1, ...  5.  Indulging  in  a  slight  abuse  of 
notation  we  can  write  the  boundary,  7,  as  an  explicit  function  of  these  5  parameters. 
More  precisely,  we  have  7  =  7u(k)U7j(k)  where  the  subscripts  u  and  l  denote  upper 
and  lower  portions  of  the  symmetric  wing  shape,  respectively  (7 i(k)  is  a  reflection 
of  7u(k)).  In  this  way  these  parameters  define  the  upper  (and  lower  by  symmetry) 
part  of  the  airfoil  by: 

7u  («)  =  {(xi,X2)\x2(xi)  =  hp(nixl  +  K2Xi  +  K3*i  +  K  4xf  +  K^?)},  (67) 

where  xlt  x2  are  the  coordinates  of  the  part  of  the  boundary  of  the  airfoil,  that  part 
lying  above  the  portion  of  the  wing  marked  by  the  c  in  the  figure  below  (see  Fig.  16) . 
Here  the  parameter  p  is  a  thickness  parameter  (for  the  NACA0012  we  take  p  —  .12). 
The  airfoil  is  extended  to  be  of  length  c'  (see  Fig.  16)  by  a  closed  trailing  edge  (see 
[24]  for  further  details  on  this  and  other  transformations,  parameterizations  and 
representations  of  the  boundary). 


C’ 

Fig.  16. 


The  NACA0012  is  described  by  the  parameter  values, 

(68)  k*  =  (0.2969,  -0.126,  -0.3516, 0.2843,  -0.1015). 

Given  parameter  values  n,  we  can  employ  the  ideas  and  methodology  from  the 
preceding  sections  of  this  paper  to  simulate  Stokes  and  Navier  Stokes  flow  over 
the  associated  k  defined  airfoil.  This  in  turn  yields  a  velocity  field  and  a  pressure 
distribution.  Let  p*  denote  the  pressure  distribution  associated  with  the  NACA0012 
values  for  n*  in  (68),  hence  p*  is  our  target  pressure  value.  Now  let  p(n)  represent 
the  pressure  distribution  and  u(k)  represent  the  velocity  field  associated  with  an 
arbitrary  set  of  parameters  k.  Similarly  we  can  write  the  airfoil  itself,  u>(k),  as  the 
region  defined  by  the  shape  associated  with  arbitrary  parameters,  k.  We  formulate 
the  inverse  problem  of  choosing  k  to  fit  p*  in,  for  example,  the  least  squares  sense, 
subject  to  certain  physical  constraints.  This  inverse  problem  can  be  stated  as  a 
general  nonlinear  programming  problem  (NLP),  namely 

min  ||p  —  P*[|z,2(7),  subject  to  K(k)  <  0.  (69) 

The  constraint  function  K(k)  in  (69)  implements  the  following  constraints: 

First,  we  require  that  (c^  -  £2|7)  <  0  for  any  x2  not  an  endpoint  of  the  airfoil. 
Here  ew  is  a  positive  minimum  airfoil  thickness  parameter  (for  instance,  we  found 


that  ew  —  10~6  works  well).  This  ensures  that  the  upper  and  lower  surface  inter¬ 
sect  in  no  more  than  two  places  (specifically  at  the  two  endpoints  of  the  airfoil). 
Secondly,  we  require  obvious  physical  upper  and  lower  bounds  on  the  volume  of  the 
airfoil,  cj(k),  the  arclength,  7,  and  the  slopes,  and 

In  practice,  the  overall  effect  of  these  constraints  is  that  they  trim  computational 
work  by  keeping  the  minimization  procedure  away  from  obviously  impossible  airfoil 
shapes  that  yield  very  small  but  non-optimal  objective  function  values.  The  im¬ 
portance  of  these  constraints  grows  if  a  large  number  of  design  variables  are  used 
in  conjunction  with  an  initial  guess  for  the  design  variables  that  results  in  an  initial 
shape  that  is  not  close  to  the  optimal  shape.  For  the  numerical  results  presented 
here,  the  constraints  resulted  in  fewer  simulations  and  hence  shorter  run-times,  but 
otherwise  did  not  drastically  affect  the  minimization  process.  This  was  in  part  due 
to  the  small  number  of  design  variables  and  the  “semi-local”  starting  value  for  the 
design  variables  (see  Fig.  17  for  a  picture  of  the  initial  and  optimal  shapes).  Fi¬ 
nally,  it  is  worth  noting  that  for  this  particular  test  problem  there  were  no  binding 
constraints  at  the  solution  (i.e.  K(k*)  <  0). 


Fig.  17:  The  starting  and  optimal  shapes  for  our  test  problems 


6.1  Optimization  Technique 

While  it  may  be  possible  to  calculate  derivatives  of  the  objective  function  in 
(69)  such  a  calculation  would  be  extremely  expensive  and  possibly  inaccurate.  An 
additional  difficulty  with  using  derivatives  in  the  solution  of  (69)  is  the  highly  non¬ 
linear  nature  of  the  objective  function.  It  can  have  many  local  minimizers  and 
inflection  points.  Our  numerical  tests  employing  derivative  based  algorithms  for 
the  solution  of  (69)  suggest  that  rapid  convergence  enjoyed  by  some  gradient  based 
methods  takes  place  only  in  an  extremely  small  neighborhood  of  the  solution,  if  at 
all.  Finally  we  would  enjoy  the  luxury  of  using  an  algorithm  that  allowed  us  to 
replace  the  smooth  least  squares  objective  function  in  (69)  with  a  continuous  but 
non-differentiable  one  (for  instance  an  L°°  objective  function). 

Algorithms  for  unconstrained  minimization  that  require  no  derivative  informa¬ 
tion  (usually  referred  to  as  direct  search  or  pattern  search  methods)  are  far  from 


new  (see,  for  example,  [27]  and  the  references  therein).  In  [28]  an  algorithm  for 
unconstrained  optimization  that  requires  no  derivative  information  was  suggested. 
A  parallel  implementation  of  the  algorithm  has  also  been  developed  and  tested  (see 
[29]).  Recently,  this  algorithm  and  its  implementation  have  been  modified  to  han¬ 
dle  constraints  (see  [30]).  The  method  samples  points  on  the  nodes  of  an  evolving 
pattern,  moving  to  the  node  in  the  pattern  with  the  value  closest  to  optimality. 
Depending  on  where  in  the  pattern  the  point  closest  to  optimality  was  located  the 
pattern  then  changes,  either  expanding  or  contracting.  When  the  pattern  expands 
its  overall  size  grows  causing  the  algorithm  to  sample  points  farther  from  its  cur¬ 
rent  location.  When  the  pattern  contracts  the  overall  size  of  the  pattern  becomes 
smaller;  this  procedure  terminates  when  the  lengths  of  the  edges  of  the  pattern  fall 
below  a  user  prescribed  tolerance. 

Pattern  search  methods  of  this  sort  typically  do  not  demonstrate  rapid  local 
convergence,  but  they  are  extremely  robust  and  far  less  susceptible  than  faster 
higher-order  methods  to  the  difficulties  introduced  when  functions  are  non-smooth 
or  the  data  is  noisy.  Pattern  search  methods  are  usually  too  slow  to  solve  optimiza¬ 
tion  problems  with  large  numbers  of  parameters.  Our  problem  (69)  is  interesting 
in  that  one  can  accurately  describe  the  shape  of  an  airfoil  with  a  relatively  small 
number  of  parameters.  The  compatibility  of  constrained  pattern  search  method 
with  our  problem  (69)  is  also  interesting.  After  numerous  numerical  tests,  we  have 
not  failed  to  locate  the  global  minimum  of  (69). 


6.2  Numerical  Example 
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Table  1:  Coarse  Minimization:  10  points  per  processor 
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Table  2:  Fine  Minimization:  25  points  per  processor 


Our  numerical  example  for  shape  optimization  requires  integrating  equations  of 
the  form  (22)  -  (23)  to  evaluate  the  objective  function.  We  consider  a  zero  degree 
and  five  degrees  angles  of  attack  for  these  test  problems  and  we  keep  the  same 
boundary  conditions  defined  by  (25)  and  (26).  Parameter  values  are  provided  in 
the  tables  below  (see  Table  1  and  Table  2).  The  cost  function  column  displays,  in 
fact,  || pc  -  P*\\l2(i)I\\p*\\l2(1)  where  pc  is  the  computed  pressure. 

We  used  16  processors  of  the  Touchstone  Delta  machine  to  perform  the  minimiza¬ 
tion  in  parallel.  A  two  stage  coarse/fine  strategy  was  employed  to  solve  the  shape 
optimization  minimization;  10  points  per  processor  were  sampled  with  a  very  large 
search  strategy  initially  in  order  to  locate  a  vague  neighborhood  of  the  solution  and 
then  25  points  per  processor  were  sampled  with  a  smaller  pattern  search  strategy 
to  polish  the  quality  of  the  solution.  For  more  details  on  constrained  pattern  search 
methods  see  [30], 

7.  CONCLUSION 

Compared  to  the  previous  results  in  [31]  it  is  clear  that  the  fictitious  domain 
methodology  that  we  advocate  has  been  substantially  progressing  and  looks  also 
promising  for  the  simulation  of  time  dependent  solution  of  viscous  flow  problems 
around  moving  rigid  bodies.  This  solution  technique  appears  to  be  well  suited 
to  the  curved  domains  on  which  we  prescribe  Dirichlet  boundary  conditions  that 
arise  from  this  class  of  shape  optimization  problems  (see  [32]  for  some  theoretical 
results).  As  always,  there  is  still  room  for  improvement.  We  mention  particularly 
the  speed  up  of  the  various  iterative  methods  used  for  the  solution  the  subproblems 
obtained  from  the  time  splitting.  Parallelization  is  also  an  important  and  timely 
issue  currently  addressed  (see  [13,  14]). 

It  appears  that  the  fictitious  domain  method  provides  an  effective  way  of  evaluat¬ 
ing  the  objective  function  for  the  Stokes  and  Navier-Stokes  flow  shape  optimization 
problem.  In  particular,  the  computational  cost  of  regriding  a  mesh  for  every  trial 
airfoil  shape  could  be  a  prohibitive  cost.  Current  work  includes  extending  these 
ideas  to  more  complicated  objective  functions  and  to  the  shape  optimization  prob¬ 
lems  associated  with  drag  reduction. 
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